Observational limit on gravitational waves from binary neutron stars in the Galaxy 
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Using optimal matched filtering, we search 25 hours of data from the LIGO 40-meter prototype 
laser interferometric gravitational-wave detector for gravitational-wave chirps emitted by coalescing 
binary systems within our Galaxy. This is the first test of this filtering technique on real interfero- 
metric data. An upper limit on the rate R of neutron star binary inspirals in our Galaxy is obtained: 
with 90% confidence, R < 0.5/hour. Similar experiments with LIGO interferometers will provide 
constraints on the population of tight binary neutron star systems in the Universe. 
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A world-wide effort is underway to test a funda- 
mental prediction of physics (the existence of gravita- 
tional waves) using a new generation of gravitational- 
wave detectors capable of making astrophysical obser- 
vations. These efforts include the US Laser Inter- 
ferometer Gravitational-wave Observatory (LIGO) |Q, 
VIRGO (French/Italian) g, GEO-600 (British/German) 
@, TAMA (Japanese) §, and ACIGA (Australian) @. 
The detectors are laser interferometers with a beam split- 
ter and mirrors suspended on wires. A gravitational wave 
displaces the mirrors, and shifts the relative optical phase 
in two perpendicular paths. This causes a shift in the in- 
terference pattern at the beam splitter Within the 
next decade, these facilities should be sensitive enough to 
observe gravitational waves from astrophysical sources at 
distances of tens to hundreds of megaparsecs (Mpc) . 

During the past 15 years, the LIGO project has used 
a 40-meter prototype interferometer at Caltech to de- 
velop optical and control elements for the full scale detec- 
tors under construction in Hanford WA and Livingston 
LA ||. In 1994, this instrument was configured as a 
modulated Fabry-Perot interferometer: light returning 
from the two arms was independently sensed |g| . In this 
configuration, the detector had its best differential dis- 
placement sensitivity of « 3.5 x 10 _19 mHz -1 / 2 over a 
bandwidth of approximately a kHz centered at 600 Hz. 

A week-long test run of the instrument was made in 
November 1994 prior to a major reconfiguration. Fig. |l| 
shows the data-taking periods. The run yielded 44.8 
hours of tape; both arms were in optical resonance for 
39.9 hours (89% of the time). Although the data was 
taken for diagnostic purposes, it provides an excellent op- 
portunity to obtain observational limits on gravitational- 
wave sources, and to examine analysis techniques. 

A major challenge arises because the real detector noise 



does not satisfy the usual simplifying assumptions: sta- 
tionary and Gaussian. The 40-m data have the expected 
colored broad-band background but with significant de- 
terministic components (spectral peaks), including ~ 10 2 
sinusoidal components arising from vibration of the sup- 
port wires and 60 Hz line harmonics. There are also 
transient features occurring every few minutes: bursts 
of noise with durations of ~ 1 — 500 ms from accidental 
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FIG. 1. Top: Boxes show data collection times. Dark 
bars show data actually filtered. Bottom: Effective distance 
D [Eq. (1)] to 90% (50%) of sources varies as the detector an- 
tenna pattern sweeps past the Galactic center. Dip at 6 hours 
is when nadir of the detector (turning with the Earth) points 
closest to the Galactic center where the potential sources are 
clustered. Fortuitously, much of the data was taken near 
such times. Jagged Line: Effective distance D at which 
a 2 x IAMq optimally oriented coalescing system would give 
SNR = p = 10. This depends on the average sensitivity of 
instrument. The small fluctuations indicate stable sensitivity. 



(natural or man-made) disturbances. These difficulties 
led us to develop data analysis techniques that make 
matched filtering methods perform well on real data. 

This Letter reports on a search of these data for binary 
inspiral chirps — the gravitational waveforms produced by 
pairs of orbiting stars or black holes. The search focuses 
on neutron star binaries in our Galaxy. On time scales 
of w 10 7 years a binary loses energy by emitting gravi- 
tational waves (primarily at twice the orbital frequency) . 
As the orbit shrinks, it circularizes and the period de- 
creases. We search for the gravitational waves that would 
be emitted during the final few seconds of this process; 
the stars orbit hundreds of times per second at separa- 
tions of tens of km before plunging together. (See || for 
results of preliminary searches.) 

The data stream was searched using matched filtering. 
This method jjj uses linear filters constructed from the 
expected waveforms, computed using the second post- 
Newtonian approximation (2PN) JsJ . The 2PN waveform 
for a 2 x 1.4M Q binary is a sweeping sinusoid which enters 
the detector pass-band around 120 Hz. The frequency 
and amplitude increase during the ensuing 255 cycles; 
after 1.35 seconds the frequency has increased to 1822 
Hz and the waveform is cut off when the stars merge. 
The 2PN approximation results in a reduction of signal 
to noise ratio (SNR) < 10% §. 

The dimcnsionless strain h(t) of the gravitational wave 
produces a differential change AL(t) = Lh(t) in the 
lengths of the two perpendicular interferometer arms [Q , 
where L = 38.25m is the average arm length. For 
a binary system (circular orbits, no spin) with masses 
M = (7774,7712) this strain is: 

h(t)= 1 PC [sina hf{t - < )+cosa hf (t - t j] . 

(1) 

Here a is a constant determined by the orbital phase and 
orientation of the binary system, to is the laboratory time 
when the chirp signal first enters the detector pass-band, 
and h^ c {t — to) are the two polarizations of the gravita- 
tional waveform produced by an inspiraling binary sys- 
tem that is optimally oriented at 1 Mpc. If z,?/-axes are 
defined by the two interferometer arms then an optimally 
oriented binary system is located on the z-axis with its 
orbital plane parallel to the x-y plane. The effective dis- 
tance D depends on the distance to the source and on 
its orientation with respect to the detector. The detec- 
tor has a non-uniform response over the sky due to its 
quadrupolar antenna pattern. If the source is not opti- 
mally oriented (i.e., not on the z-axis or the orbital plane 
is tipped), then D is greater than the source-detector dis- 
tance. The formulae for h^f c are Eqs. (2, 3a, 4a) of Ref. ||. 

The detector signal is the voltage applied to produce 
a feedback force on the mirrors to hold the interferom- 
eter in resonance; it is proportional to the differential- 
displacement AL(t). This voltage v(t) was recorded at 



a sample rate of 9868.42 Hz by a 12 bit analog-to-digital 
converter. Quantizing the data reduces the SNR by less 
than 0.9% pj. The instrument's frequency and phase 
response R(f) was determined at the beginning of each 
of eleven ~ 4 hour data runs by applying known pertur- 
bative forces to the interferometer jO]. These eleven 
calibration curves differ by less than 5%. Because errors 
in calibration affect the SNR only at second order, we 
estimate the effects of any calibration errors or drifts on 
SNR to be less than 0.3%. The voltage output Vh(t) that 
would be produced by a binary inspiral is given by 

v h (t)= I R(t-t')h(t')dt' 

J — OO 

/OO 
-00 

where Q(t) and Q(f) denote Fourier-transform pairs. 

We search for inspiral waveforms using (digital) 
matched filtering. Because the inspiral waveforms de- 
pend upon the source masses M = (7711,7772) we use 
a "bank" of template waveforms with masses spaced 
closely enough [jl2| to detect any signal in the mass range 
1.0 M Q < mi, m 2 < 3.0 M© 0]. The bank contains 687 
filters Mfc and is designed so that no more than 2% of 
SNR would be lost if the mass parameters M of a signal 
did not exactly match one of the Mk . For each mass pair 
Mfc in the template bank two real signals are constructed: 

Y s,c ( .\ _ AT s,c f v(f)h s , c k R(f) 2wift jr (o\ 

Xk {t) - N " U s v {\f\) e df - (2) 

These are the outputs of optimal filters matched to the 
waveform of the fcth mass-pair The denominator 

Sv{\f\) is (an estimate of) the one-sided power spectral 
density of v(t); if the detector's noise is stationary and 
Gaussian, then these filters are optimal. The normaliza- 
tion factor N^' c is chosen so that, in the absence of any 
signals, the mean value of [X k ' c (t)} 2 is unity. We define 
the SNR for the kth template waveform to be 

Pk (t) = SNR = ^[X> k {t)f + [X<(t)] 2 , 

arrived at by maximizing over the phase a of the binary 
system. The effective distance D at which coalescence 
of 2 x 1.4M Q stars would yield an SNR of 10 in the in- 
terferometer is shown in Fig. El (The definition of the 
SNR follows Ref. |E| and other literature. Its expected 
value for a source scales oc D" 1 . Its rms value for a single 
template is y2 in the presence of Gaussian noise alone.) 

The data was processed, using FFT methods, in over- 
lapping w 26.6s segments (2 18 samples). To avoid end 
effects, S^d/I) in Eq. (g) was truncated at w 13.3s in 
the time domain. The longest chirp signal was w 2.4 s 
long, so the data were overlapped by the total filter im- 
pulse response time of ss 15.6 s (155 072 samples) giv- 
ing a filter output duration of « 10.85 s/segment. Since 



the process of bringing the optical cavities into resonance 
(lock) excites vibrations in the suspension wires, we dis- 
carded the first three minutes of data after each lock 
acquisition, allowing the vibrations to damp below other 
noise sources. Of the 39.9 locked hours of data, 8.8 hours 
were in intervals too short to analyze; 111 locked inter- 
vals remained. Discarding the startup transient impulse 
response of the filters and the first three minutes of lock 
yielded 39.9 — 8.8 — 6.0 = 25.0 hours of data analyzed, in 
8289 intervals of filter output (top of Fig. ||) . 

Poorly understood, non-stationary noise events cor- 
rupt the data. However, these transient events do not 
have the time-frequency behavior of inspiral chirps, so 
we can use the broad-band nature of the interferomet- 
ric detector to reject them. These events are discrimi- 
nated from chirps by a x 2 time-frequency test (Sec. 5.24 
of Ref. 13 1). The frequency band (DC to Nyquist) is 
divided into p subintervals, chosen so that for a chirp 
superposed on Gaussian noise with the observed power 
spectrum the expected contribution to p is equal for each 
subinterval. One forms a statistic x 2 by summing the 
squares of the deviations of the p signal values from the 
expected value for the two template polarizations. We 
choose p — 20 so that Galactic signals that fall at the 
maximal template mismatch would not be rejected. In 
the presence of Gaussian noise plus chirp the statistic has 
a x 2 distribution with 2p—2 = 38 degrees of freedom [ Oj . 

Occasionally, there are short sections (i.e., glitches) 
in the data when the instrument's output significantly 
exceeds the rms value. Some of these glitches were 
seismically-induced. These short sections cause the out- 
puts of the optimal filters to ring, but do not resem- 
ble binary inspiral chirps and are uniformly rejected by 
the time-frequency technique described above. However 
these glitches bias S v (\f\) enough to create non-optimal 
filters. To prevent this problem we estimate the power 
spectrum by averaging it for the 8 glitch-free segments 
closest in time to the section being analyzed. The glitches 
were identified by seeing if too many samples fell outside 
a ±3cr range or any fell outside a ±5<7 range. The num- 
ber of segments (8) was chosen to reduce the variance of 
the spectrum while still tracking changes in instrument 
performance. 

The data was processed in about 32 hours of clock 
time on a 48 node Beowulf computer at UWM (29 Gflops 
peak). The output of the filtering process is a list of 
signals for each segment j: the maximum (over t) SNR 
obtained for each filter k in the bank of 687 filters, the 
time tj at which that maximum occurred, the value of the 
X 2 statistic for that filter, and N s,c . In a given segment of 
data, we say that an event has occurred if the maximum 
SNR, over all filters for which the statistic x 2 ues below 
some threshold x*? exceeds a threshold p*. The total 
number N of events observed in the data set of Fig. |l| is 
plotted as a function of these thresholds in Fig. ||[ 

Without operating two or more detectors in coinci- 
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FIG. 2. Top: total number N of events observed, as 
a function of the SNR threshold p„ and the threshold x*- 
Bottom: fraction e of Galactic inspiral chirp signals that 
would lie above SNR threshold p* and below y- 2 threshold x* ■ 



dence, it is impossible to characterize the non-Gaussian 
and non-stationary background well enough to state with 
confidence that an event has been detected. However one 
may estimate upper limits on the rate of Galactic neutron 
star binary inspirals (Poisson-distributcd in time) using 
a method which requires minimal assumptions about de- 
tector noise. Our limit i?go% is based on the probability 
of a Galactic neutron star binary signal having an SNR 
as big as the largest SNR observed. If the actual inspiral 
rate is greater than i?9o% then it is likely that we would 
have observed a larger SNR event. Fig. |l| shows that 
much of the time the detector was not pointing at the 
Galactic bulge; therefore the detector was only sensitive 
to a fraction of Galactic binary inspirals. Thus the event- 
rate bound depends on two numbers: (i) the efficiency 
£max with which the instrument and filtering/analysis 
process can detect a binary inspiral in the Galaxy at the 
SNR pmax of the largest observed event, and (ii) the total 
length T = 25.0 hours of filtered data. 

We determined the efficiency e by Monte-Carlo simu- 
lation, doing additional runs through the data set, and 
adding simulated Galactic inspiral waveforms [convolved 
with the detector response function R(f)] at 30 s intervals 
into the detector output v(t). This allows us to charac- 
terize the detection process with the properties of the real 
instrument noise rather than an ad-hoc model. The in- 
serted waveforms were drawn from a population of binary 
neutron stars with a spatial number distribution given by 
dN oc e~ v2 / 2V «VdV x e -\ z \/ ,lz dZ where V is Galacto- 
centric radius, Vq = 4.8 kpc, Z is height off the Galactic 
plane, and hz = 1 kpc is the scale height. This distri- 
bution is similar to the one presented in Ref. 14|. 



The 

detection efficiency e is the fraction of these simulated in- 
spirals which registered as events in our filtering/analysis 
procedure; it increases as the SNR threshold is de- 
creased, or as x* is increased, and is shown in Fig. || for 



the most-probable mass range Jll| of 1.29 to 1.45 M Q 
(the results depend weakly on the mass). 

Our analysis gives an event rate bound. With 90% 
confidence, the rate of binary inspirals in our Galaxy is 
less than #90% = 3.89/[Te(/5 max , x*)] where p max = 8.34 
is the largest SNR event observed, and the threshold 
xl = 49.5 is chosen so that there is a 10% chance of 
rejecting a real chirp signal in stationary Gaussian noise. 
(This is a Bayesian credible interval computed using a 
"uniform prior" for the event rate. The dimensionlcss 
numerator depends only on the confidence level.) The 
efficiency e(8.34, 49.5) = 0.33 gives Rqq% — 0.5/hour. 
This is a 90% confidence limit if the largest event is a 
real binary inspiral event. If the largest event is noise, 
the confidence is > 90%. Thus, i?go% gives a conservative 
upper limit on the event rate when the detector noise is 
poorly understood. [The on-site environmental monitors 
show that some of the larger events in Fig. || arise from 
seismic disturbances or laser power fluctuations, but the 
largest event (on which our rate limit is based) was de- 
tected during normal instrument operation.] 

Let us compare our limit i?go% = 0.5/hour, with the 
limit that could be obtained from the ideal analysis of 
an instrument that could detect every Galactic event. 
Operating for the same total time T = 25.0 hours with 
an efficiency e = 1, the limit obtained would be three 
times better: i?go% = 0.17/hour. 

Using stellar population models |l6[ ], one can forecast 
an expected inspiral rate of R ~ lO^^r -1 , far below our 
limit. However, unlike these model-based forecasts, our 
inspiral limit is based on direct observations of inspirals. 
Our study also demonstrates methods being developed 
to analyze data from the next generation of instruments. 

A previous search using 100 hours of coincident Glas- 
gow/Garching interferometer data gave an upper limit on 
burst sources jl7). The current generation of resonant- 
mass detectors [|l8| has established upper limits on 
monochromatic signals and stochastic background, but 
neither search addressed the binary inspiral rate. A coin- 
cidence analysis of bar data for coalescing binaries might 
produce a stronger limit than ours. 

The full-scale 4- km LIGO interferometers will be much 
more sensitive than the 40-meter prototype. Comprehen- 
sive instrument monitoring will permit detailed charac- 
terization of instrument anomalies and removal of some 
environmental noise. Correlation between three indepen- 
dent instruments will provide lower false alarm rates and 
greater statistical confidence. This will augment the tech- 
niques used here and allow LIGO to detect sources, as 
well as set tight rate limits. For example, if the largest 
coincident event detected by the LIGO interferometers 
has a SNR p max = 5.5, then we would obtain the limit 

^o^exio-Mpc-V" 1 ^) 3 ^ , 

\ ''max / \ ^ obs / 

on the rate of inspiral in the universe, where T b s is the 



observation time, and r mm is the distance to an opti- 
mally oriented source with SNR /O max = 5.5. For the ini- 
tial LIGO interferometers, the distance is r max = 55 Mpc; 
it will be ten times larger for the enhanced interferome- 
ters, giving an expected rate limit of 6 x 10 _8 Mpc _3 yr _1 . 
These limits should be compared to the best guess rate 
of 8 x 10~ 8 Mpc _3 yr _1 given by Phinney Q. 
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